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In spite of over a century of research on cortical circuits, it is still unknown how 
many classes of cortical neurons exist. In fact, neuronal classification is a difficult 
problem because it is unclear how to designate a neuronal cell class and what are the 
best characteristics to define them. Recently, unsupervised classifications using cluster 
analysis based on morphological, physiological, or molecular characteristics, have provided 
quantitative and unbiased identification of distinct neuronal subtypes, when applied to 
selected datasets. However, better and more robust classification methods are needed 
for increasingly complex and larger datasets. Here, we explored the use of affinity 
propagation, a recently developed unsupervised classification algorithm imported from 
machine learning, which gives a representative example or exemplar for each cluster. 
As a case study, we applied affinity propagation to a test dataset of 337 interneurons 
belonging to four subtypes, previously identified based on morphological and physiological 
characteristics. We found that affinity propagation correctly classified most of the neurons 
in a blind, non-supervised manner. Affinity propagation outperformed Ward's method, a 
current standard clustering approach, in classifying the neurons into 4 subtypes. Affinity 
propagation could therefore be used in future studies to validly classify neurons, as a first 
step to help reverse engineer neural circuits. 



Keywords: affinity propagation, cortex, interneurons, cell types 



INTRODUCTION 

To properly understand the structure and function of any circuit 
it seems essential to objectively define its elements. Unfortunately, 
as opposed to elements in electronic circuits, neurons in brain 
circuits do not come pre-labeled and it is not clear exactly what 
comprises a neuronal cell type. GABAergic neocortical interneu- 
rons are a particularly difficult case, due to their large molecular, 
morphological and physiological diversity (Fairen et al., 1984; 
Mott and Dingledine, 2003; Ascoli et al, 2008). In the past, cell 
type classification was a qualitative and subjective task that led 
to inconsistent classes of neurons. Recently, quantitative meth- 
ods using unsupervised cluster analysis have become standard 
for classification of neurons (Cauli et al, 1997; Karube et al., 
2004; Ma et al, 2006; Dumitriu et al., 2007; Helmstaedter et al., 
2009; Karagiannis et al, 2009; McGarry et al, 2010; DeFelipe 
et al., 2013). In particular, traditional cluster analysis using Ward's 
method has been effective, but it is a simple technique with some 
drawbacks. Hierarchical agglomerative clustering is a bottom-up 
technique, that is, it starts by grouping the two "closest" cells as 
defined by the algorithm, then joins the next "closest" sets and 
so on. As a hierarchical clustering technique the grouping deci- 
sions it makes are inflexible so, once two cells are joined together 
they remain joined in the final hierarchy. Moreover, hierarchical 
methods are susceptible to a chaining effect, where cells are some- 
times assigned to existing clusters rather than being grouped in 
new clusters. These qualities of the method could pose limitations 



as it prevents it from testing multiple possible groupings of the 
dataset. 

A new, and more sophisticated, classification method called 
affinity propagation does not have those limitations, is an unsu- 
pervised algorithm, and has demonstrated success with greatly 
improved results over standard methods in other classification 
problems (Frey and Dueck, 2007). In particular, affinity prop- 
agation goes beyond the solution of the neuronal classification 
problem and solves a related but also difficult problem: the iden- 
tification of exemplars, i.e., elements that best characterize each 
class (Mezard, 2007). In neuroscience, affinity propagation was 
recently used for the analysis of spike synchronization in rat 
brains (Takahashia et al., 2010), in learning the role of sleep slow 
wave activity in visuomotor learning (Landsness et al., 2009), and 
in the discovery of synaptic connectivity organizational princi- 
ples (Perin et al., 2011). When applied to neuronal classifications, 
affinity propagation could therefore identify exemplar neurons 
to be used as a compressed, characteristic representation of each 
neuron subtypes. This seems advantageous, since focusing the 
analysis on exemplars could ease the burden of analyzing hun- 
dreds, or thousands of neurons, for the identification of salient 
characteristics of neuron classes. 

In this work we explored the application of affinity prop- 
agation to neuronal classification, by using the algorithm to 
blindly classify a test dataset of four known interneuron sub- 
types. The test dataset was comprised of 67 morphological 
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FIGURE 1 | Morphological Clusters. Neurons are represented as colored neuron type. BC and ChC (red and magenta) are closely related PV+ 

glyphs. Colors red, blue, green, and magenta, respectively, represent neuron interneurons and MC and non-MC (blue and green) are both subtypes of 

types BC, MC, non-MC, and ChC. The label of the exemplar in each cluster is SOM+ cells. Note how cluster 3 (from top) groups BC and ChC jointly and 

shaded in yellow. Ten clusters are found; most clusters are dominated by a cluster 10 (last one) groups MC and non-MC together. 



Table 1 | Accuracies computed for the Morphology, Physiology, and 


Morphology + Physiology databases. 




Morphology 


Physiology 


Morphology-!- Physiology 


database 


database 


database 


ncluster ACC 


ncluster ACC 


ncluster ACC 


10 0.7374 


36 0.8505 


8 0.7857 



ncluster is the number of clusters found by the algorithm. ACC, Accuracy 
obtained using the four known classes of neurons (BC, ChC, MC, non-MC). 



and 20 electrophysiological variables (Supplementary Table 1) 

describing (1) parvalbumin-positive (PV+) basket cells (BC), (2) 
PV+ chandelier cells (ChC), (3) somatostatin-positive (SOM+) 
Martinotti cells (MC), and (4) SOM+ non-Martinotti cells 
(non-MC), as described in (McGarry et al., 2010). We found 
that affinity propagation generates an accurate classification in 
separating these four known interneuron subtypes. Our data 
suggest that affinity propagation could be a powerful new 
classification tool for discovering or defining neuronal cell 
types. 

MATERIALS, METHODS, AND PROTOCOL 
PREPARATION OF BRAIN SLICES 

Acute brain slices were prepared from Nkx2.1, G42, or GIN 
mice, with an average age of 15 postnatal days (range P13- 
25). Mice were quickly decapitated, the brain was removed 
and then immediately placed in cold sucrose cutting solu- 
tion (222 mM sucrose, 2.6 mM KC1, 27 mM NaHC0 3 , 1.5 mM 
NaH 2 P0 4 , 0.5 mM CaCl 2 , 3 mM MgSQ 4 , bubbled with 95% 0 2 , 



5%C02). Coronal slices of 300 [im thickness were cut using a 
Vibratome and transferred to a holding chamber at room tem- 
perature with oxygenated ACSF (126 mM NaCl, 3 mM KC1, 3 mM 
MgS0 4 , 1 mM CaCl 2 , 1.1 mM NaH 2 P0 4 , 26 mM NaHC0 3 , and 
10 mM dextrose, bubbled with 95% 0 2 , 5%C0 2 ). After at least 
an additional 30min equilibration at room temperature, slices 
were transferred to a recording chamber with perfusion of ACSF 
bubbled with 95% 0 2 , 5%C0 2 . 

TRANSGENIC MOUSE LINES 

To identify different types of interneurons we used three trans- 
genic mouse lines. First, we used the G42 line that labels PV+ 
cells (Chattopadhyaya et al, 2007). PV+ cells are fast spiking 
interneurons with basket or ChC morphology. We could iden- 
tify BC from ChC by their distinctive morphologies and threshold 
spiking responses. Additionally, to find ChCs we used Nkx2.1 Cre 
MADM mice (referred to as Nkx2.1 mice) where we could often 
identify ChC by their distinct axonal arbors, visible with illumi- 
nation of the GFP (Woodruff et al, 2009). The Nkx2.1 line labels 
a population of interneurons that express the transcription factor 
Nkx2.1, labeling interneurons that migrate from the medial gan- 
glionic eminence (MGE), thus including ChCs (Xu et al., 2008). 
In both G42 and Nkx 2.1 lines, we had a high success rate of find- 
ing ChC at the top of layer 2, close to the layer 1 border, where 
the normally rare ChC were found among GFP+ cells with a 
probability of 50-70% (Woodruff et al., 2009). Finally, we used 
the GIN line to label SOM+ cells (Oliva et al, 2000). SOM+ 
cells are regular spiking interneurons with varied morphology. 
We previously identified three subtypes of SOM+ interneurons 
in GIN mice based on morphology and physiology-MC and 
two novel subtypes (McGarry et al, 2010). Following that work, 
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FIGURE 2 | Physiological Clusters. Code as in Figure 1. See text for details on the clusters. 
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FIGURE 3 | Clusters of the combined Morphology+Physiology database. Code as in Figure 1 . See text for details on the clusters. 



here we distinguish between the MC and two novel subtypes 
(non-MC). 

ELECTR0PHYSI0L0GY RECORDINGS 

Slices were placed in a recording chamber at room temperature 
with flowing oxygenated ACSF. Pipettes of 3-7 MQ, resistance 
were pulled from borosilicate glass. Whole cell recordings were 
taken in current clamp mode. Only cells with healthy resting 
membrane potential (between —55 and — 80 mV) were selected 
for recording. Supplementary Figure 2 shows a representative set 
of the traces of the complete dataset. 

ELECTR0PHYSI0L0GY ANALYSIS 

Twenty variables were measured for each neuron by analysis of 
the recordings in MATLAB (Supplementary Table 1). Variables 
describing firing and passive properties were based on the Petilla 
terminology (Ascoli et al., 2008). 

HISTOLOGICAL PROCEDURE 

Neurons were filled with biocytin by a patch pipette. Slices 
were kept overnight in 4% paraformaldehyde in 0.1 M phos- 
phate buffer (PB) at 4°C. Slices were then rinsed three times for 
5min per rinse on a shaker in 0.1 M PB. They were placed in 
30% sucrose mixture (30 g sucrose dissolved in 50 ml ddF^O and 



50 ml 0.24 M PB per 100 ml) for 2h and then frozen on dry ice 
in tissue freezing medium. The slices were kept overnight in a 
— 80°C freezer. The slices were defrosted and the tissue freezing 
medium was removed by three 20 min rinses in 0.1 M PB. Slices 
were kept in 1% hydrogen peroxide in 0.1 M PB for 30 min to 
pretreat the tissue, then were rinsed twice in 0.02 M potassium 
phosphate saline (KPBS) for 20 min. The slices were then kept 
overnight in Avidin-Biotin-Peroxidase Complex. Next the slices 
were rinsed three times in 0.02 M KPBS for 20 min each. Each 
slice was then placed in DAB (0.7 mg/ml 3,3"-diaminobenzidine, 
0.2 mg/ml urea hydrogen peroxide, 0.06 M Tris buffer in 0.02 M 
KPBS) until the slice turned light brown, then immediately trans- 
ferred to 0.02 M KPBS and transferred again to fresh 0.02 M KPBS 
after a few minutes. The stained slices were rinsed a final time in 
0.02 M KPBS for 20 min. Each slice was observed under a light 
microscope and then mounted onto a slide using crystal mount. 

RECONSTRUCTION AND ANALYSIS OF MORPHOLOGY 

Successfully filled and properly stained neurons were then 
reconstructed using Neurolucida software (MicroBrightField) 
(Supplementary Figure 1). The neurons were viewed with a 
lOOx oil objective on an Olympus 1X71 inverted light micro- 
scope or an Olympus BX51 upright light microscope. The 
neuron's processes were traced manually while the program 
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Table 2 | Significant variables for each morphological cluster. 



Variables 


CI1 


CI2 


CI3 


CI4 


CI5 


CI6 


CI7 


CI8 


CI9 


CI10 


Number of dendrites 




















Dendritic node total 




















Dendrite node density 




















Total dendritic length 




















Average length of dendrites 




















Total surface area of dendrites 
















+ 




Ratio of dendritic length to surface area 




















Dendritic torsion ratio 




















Dendritic planar angle ave 




















Dendritic planar angle stdv 




















Dendritic local angle ave 
















_ 




Dendritic local angle stdv 
















_ 




Dendritic spline angle ave 
















— 




Dendritic spline angle stdv 




















Ave tortuosity of dendritic segments 




+ 
















Stdv of tortuosity of dendritic segments 




















Dendritic segment length ave 




















Dendritic segment length stdv 




















Ave tortuosity of dendritic nodes 




















Stdv tortuosity of dendritic nodes 




















Number of dendritic sholl sections 




















Dendritic sholl length at 50 \im 




_ 














■ 


_ 


Dendritic sholl length at 100 \im 




— 


+ 














— 


Dendritic sholl length at 1 50 |j,m 






















Convex hull dendrite area 




















Convex hull dendrite perimeter 




















Convex hull dendrite volume 
















■ 




Convex hull dendrite surface area 
















+ 




Highest order dendritic segment 




















k-dim dendrites (fractal analysis) 




















Axonal node total 












_ 




+ 




Total axonal length 












_ 




+ 




Ratio of axonal length to surface area 


+ 


+ 


















Highest order axon segment 


+ 












_ 








Axonal torsion ratio 




















Axonal planar angle ave 




















Axonal planar angle stdv 




















Axonal local angle ave 
















_ 




Axonal local angle stdv 




















Axonal spline angle ave 




















Axonal spline angle stdv 




















Ave tortuosity of axonal segments 




















Stdv of tortuosity of axonal segments 




















Axonal segment length ave 














■ 








Axonal segment length stdv 














+ 








Ave tortuosity of axonal nodes 


+ 


















Stdv tortuosity of axonal nodes 


+ 


















Number axonal sholl sections 










+ 










Axonal sholl length at 100|j,m 








■ 














Axonal sholl length at 200 \i.m 
















■ 





(Continued) 
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Table 2 | Continued 



Variables 


CI1 


CI2 


CI 3 


CI4 


CIS 


CI6 


CI7 


CIS 


CI9 


CI10 


Axonal sholl length at 300 


+ 






— 














Axonal length density2 












— 




+ 




Axonal node density 


+ 












— 








Axonal node density2 












— 




+ 




Convex hull axon area 






— 




+ 


— 




+ 




Convex hull axon perimeter 






— 




+ 










Convex hull axon volume 






— 




+ 


— 




+ 




Convex hull axon surface area 






— 




+ 


— 




+ 




k-dim axon (fractal analysis) 
















+ 




Total surface area of axon 
















+ 




Somatic perimeter 




















Somatic area 














+ 






Somatic aspect ratio 




















Somatic compactness 




















Somatic form factor 




















Somatic roundness 




















Relative distance to pia 























The "+ " or "— " signs mean that the average value of a variable in a cluster is significantly higher (or lower), than the average value of that variable in the database. 



recorded the coordinates of the tracing to create a digital three- 
dimensional reconstruction. In addition to the neuron, the pia 
and white matter were drawn. The Neurolucida Explorer program 
was used to measure 67 morphological variables of the recon- 
struction describing somatic, dendritic, and axonal properties. 
(Supplementary Table 1). 

AFFINITY PROPAGATION 

Affinity propagation is a clustering algorithm that, given a set 
of points and a set of similarity values between the points, 
finds clusters of similar points, and for each cluster gives a 
representative example called an exemplar (Frey and Dueck, 
2007). Affinity propagation has several advantages over related 
techniques. Methods such as fc-centers clustering and fc-means 
clustering store a relatively small set of estimated cluster cen- 
ters at each step. These techniques can be improved by using 
methods that begin with a large number of clusters and then 
prune them, but they still rely on random sampling and 
make hard pruning decisions that cannot be recovered from. 
In contrast, by simultaneously considering all data points as 
candidate centers and gradually identifying clusters, affinity prop- 
agation is able to avoid many of the poor solutions caused 
by unlucky initializations and hard decisions (Frey and Dueck, 
2007). 

A characteristic that makes affinity propagation different from 
other clustering algorithms is that the points directly exchange 
information between them regarding the suitability of each point 
to serve as an exemplar for a subset of other points. The algo- 
rithm takes as input a matrix of similarity measures between each 
pair of points s(i, k). Instead of requiring that the number of clus- 
ters be predetermined, affinity propagation takes as input a real 
number s(k, k) for each data point k. These values, which are 
called preferences, are a measure of how likely each point is to 



be chosen as exemplar. In our case the preference can be under- 
stood as a particular weight given to each neuron according to a 
priori knowledge of the suitability of the neurons to be exemplars. 
This parameter can be used to bias the clustering procedure when 
there are some neurons that are known to be good descriptors. 
However, in our experiments we did not assume any a priori 
information and all neurons where given the same preference. 
The algorithm works by exchanging messages between the points 
until a stop condition, which reflects an agreement between all 
the points with respect to the current assignment of the exem- 
plars, is satisfied. These messages can be seen as the way the 
points share local information in the gradual determination of the 
exemplars. 

There are two types of messages to be exchanged between data 
points. The responsibility r(z, k), sent from data point i to can- 
didate exemplar point k, reflects the accumulated evidence for 
how well-suited point k is to serve as the exemplar for point i, 
taking into account other potential exemplars for point i. The 
availability a(i, k), sent from candidate exemplar point k to point 
i, reflects the accumulated evidence for how appropriate it would 
be for point i to choose point k as its exemplar, taking into 
account the support from other points that point k should be an 
exemplar. 

The availabilities are initialized to zero: a(i,k) = 0. Then, the 
responsibilities are computed using the rule: 

r(i, k) <- s(z, k) - max^ij^i; {a(i, k') + s(i, k')) (1) 

The responsibility update shown in Equation (1) allows all the 
candidate exemplars compete for ownership of a data point. 
Evidence about whether each candidate exemplar would make a 
good exemplar is obtained from the application of the following 
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Table 3 | Significant variables for each physiological cluster. 



Variables CM 


CI2 


CI3 


CI4 


CI5 


CI6 


CI7 


CI8 


CI9 


CI10 


Rheobase (pA) 




















Resting membrane potential (mV) 


■ 












■ 






AP1 Amplitude (mV) 














+ 






AP1 duration (ms) 


+ 


+ 


■ 




■ 


+ 


+ 






AP1 half-width (ms) 


+ 




+ 




+ 


+ 


+ 






AP1 rise time (ms) + 


+ 


+ 


+ 




+ 


+ 


+ 






AP1 fall time (ms) 


+ 


+ 


+ 




+ 


+ 


+ 






AP1 rise rate (mV/ms) 


— 




— 




_ 




— 






AP1 fall rate (mV/ms) 




















AP2 Amplitude (mV) 














■ 






AP2 duration (ms) 


■ 




■ 




■ 


■ 


+ 






AP2 half-width (ms) 


+ 




+ 




+ 


+ 


+ 






AP2 rise time (ms) + 


+ 


+ 


+ 




+ 


+ 


+ 






AP2 fall time (ms) 


+ 


+ 


+ 




+ 


+ 


+ 






AP2 rise rate (mV/ms) 




















AP2 fall rate (mV/ms) 




















AP Drop (mV) 




















Index of accommodation + 




















Input resistance 


■ 




■ 




■ 




■ 






Spike frequency — 




_ 
















cm 


CM 2 


CI13 


C1 14 


CI15 


C 11 6 


C1 17 


CI18 


C1I9 


CI20 


Rheobase (pA) 


+ 
















+ 


Resting membrane potential (mV) 


— 
















— 


AP1 Amplitude (mV) 








_ 




_ 


_ 


_ 




AP1 duration (ms) 




















AP1 half-width (ms) 




















AP1 rise time (ms) 




















AP1 fall time (ms) 


— 


















AP1 rise rate (mV/ms) 




















AP1 fall rate (mV/ms) 


+ 


















AP2 Amplitude (mV) 








_ 




_ 








AP2 duration (ms) 


_ 


















AP2 half-width (ms) 




















AP2 rise time (ms) 




















AP2 fall time (ms) 




















AP2 rise rate (mV/ms) 


■ 


















AP2 fall rate (mV/ms) 


+ 


















AP Drop (mV) 


+ 
















■ 


Index of accomodation 


















+ 


Input resistance 


















_ 


Spike frequency 




















CI21 


CI22 


CI23 


CI24 


CI25 


CI26 


CI27 


CI 28 


CI29 


CI30 


Rheobase (pA) 












+ 






+ 


Resting membrane potential (mV) 










+ 










AP1 Amplitude (mV) 




















AP1 duration (ms) 




















AP1 half-width (ms) 




















AP1 rise time (ms) 





















(Continued) 
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Table 3 | Continued 



Variables CI1 


CI2 


CI3 


CI4 


CI5 


CI6 


CI7 


CI8 


CI9 


CI10 


AP1 fall time (ms) 




















AP1 rise rate (mV/ms) 


+ 






+ 


+ 










AP1 fall rate (mV/ms) 




+ 




+ 


+ 










AP2 Amplitude (mV) 








+ 


+ 










AP2 duration (ms) 




















AP2 half-width (ms) 




















AP2 rise time (ms) 




















AP2 fall time (ms) 




















AP2 rise rate (mV/ms) 


+ 






+ 


+ 










AP2 fall rate (mV/ms) 




+ 




+ 


+ 










AP Drop (mV) 












+ 








Index of accomodation 




















Input resistance 




















Spike frequency 




+ 








+ 








CI31 


CI32 


CI33 


CI34 


CI35 


CI36 










Rheobase (pA) 


+ 


















Resting membrane potential (mV) 




















AP1 Amplitude (mV) 




















AP1 duration (ms) 




















AP1 half-width (ms) 




















AP1 rise time (ms) 




















AP1 fall time (ms) 




















AP1 rise rate (mV/ms) 






+ 














AP1 fall rate (mV/ms) 






+ 














AP2 Amplitude (mV) 




















AP2 duration (ms) 




















AP2 half-width (ms) 




















AP2 rise time (ms) 




















AP2 fall time (ms) 




















AP2 rise rate (mV/ms) 




















AP2 fall rate (mV/ms) 




















AP Drop (mV) 




















Index of accomodation 




















Input resistance 




















Spike frequency 





















availability update: 



a{i, k) <<-min|o, r(k, k) + ^ f . k max {o, r(/, k)} J (2) 



In the availability update shown in Equation (2) only the positive 
portions of incoming responsibilities are added, because it is only 
necessary for a good exemplar to explain some data points (posi- 
tive responsibilities), regardless of how poorly it explains points 
with negative responsibilities. To limit the influence of incom- 
ing positive responsibilities, the total sum is thresholded so that 
it cannot go above zero. 



The self-availability a(k, k) is updated differently: 

a(k, k) jt max {°' r (»', *)} ( 3 ) 

For a point i, the value of k that maximizes a(i, k) + r(i, k) either 
identifies point i as an exemplar if k = i, or identifies the data 
point that is the exemplar for point i. 

Update rules described by Equations (1), (2), and (3) require 
only local computations. The message-passing procedure may be 
terminated after a fixed number of iterations, when changes in 
the messages fall below a threshold, or after the local decisions 
stay constant for some number of iterations. 
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Table 4 | Significant variables for each cluster of the combined Morphology+Physiology database. 



Variables CM 


CI2 


CI3 


CI4 


CI5 


CI6 


CI7 


CI8 


Number of dendrites 
















Dendritic node total 
















Dendrite node density 
















Total dendritic length 
















Average length of dendrites 
















Total surface area of dendrites 
















Ratio of dendritic length to surface area 
















Dendritic torsion ratio 
















Dendritic planar angle ave 
















Dendritic planar angle stdv 
















Dendritic local angle ave 
















Dendritic local angle stdv 
















Dendritic spline angle ave 
















Dendritic spline angle stdv 
















Ave tortuosity of dendritic segments 
















Stdv of tortuosity of dendritic segments 
















Dendritic segment length ave 
















Dendritic segment length stdv 


_ 














Ave tortuosity of dendritic nodes 
















Stdv tortuosity of dendritic nodes 
















Number of dendritic sholl sections 
















Dendritic sholl length at 50 \im 








+ 




_ 




Dendritic sholl length at 100 |xm 








+ 








Dendritic sholl length at 1 50 |j,m 
















Convex hull dendrite area 
















Convex hull dendrite perimeter 




+ 












Convex hull dendrite volume 
















Convex hull dendrite surface area 
















Highest order dendritic segment 
















k-dim dendrites (fractal analysis) 
















Axonal node total 












■ 




Total axonal length 












+ 




Ratio of axonal length to surface area 
















Highest order axon segment 












+ 




Axonal torsion ratio 
















Axonal planar angle ave 
















Axonal planar angle stdv 
















Axonal local angle ave 
















Axonal local angle stdv 
















Axonal spline angle ave 
















Axonal spline angle stdv 
















Ave tortuosity of axonal segments 
















Stdv of tortuosity of axonal segments 
















Axonal segment length ave 




+ 








— 




Axonal segment length stdv 




+ 












Ave tortuosity of axonal nodes 
















Stdv tortuosity of axonal nodes 
















Number axonal sholl sections 
















Axonal sholl length at 100|j,m 
















Axonal sholl length at 200 |j,m 
















Axonal sholl length at 300 i^m 
















Axonal length density2 












■ 





(Continued) 
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Table 4 | Continued 



Variables CM 


CI2 


CI3 


CI4 


CI5 


CI6 


CI7 


CI8 


Axonal node density 












+ 




Axonal node density2 












+ 




Convex hull axon area 
















Convex hull axon perimeter 
















Convex hull axon volume 
















Convex hull axon surface area 
















k-dim axon (fractal analysis) 












+ 




Total surface area of axon 
















Somatic perimeter 
















Somatic area 
















Somatic aspect ratio 
















Somatic compactness 
















Somatic form factor 
















Somatic roundness 
















Relative distance to pia 












_ 




Rheobase (pA) 








— 








Resting Membrane Potential (mV) 
















AP1 Amplitude (mV) 
















AP1 duration (ms) 








+ 




_ 




AP1 half-width (ms) 








+ 




_ 




AP1 rise time (ms) 








+ 




— 




AP1 fall time (ms) 








+ 








AP1 rise rate (mV/ms) 












■ 




AP1 fall rate (mV/ms) 












+ 




AP2 Amplitude (mV) 
















AP2 duration (ms) 








■ 




— 




AP2 half-width (ms) 








+ 




_ 




AP2 rise time (ms) 








+ 








AP2 fall time (ms) 








+ 








AP2 rise rate (mV/ms) 












■ 




AP2 fall rate (mV/ms) 












+ 




AP Drop (mV) 
















Index of accomodation 
















Input resistance 
















Spike frequency 

















Similarly to other propagation methods, damping should be 
used to confront numerical oscillations that arise in some cir- 
cumstances. This technique consists of setting each message to 
X times its value from the previous iteration plus 1 — "k times 
its prescribed updated value (0 < X < 1). A pseudocode of our 
approach for neuron classification using the affinity propagation 
algorithm is shown in Algorithm 1. 



Algorithm 1 : Neuron classification using affinity propagation 

1 . Normalize each of the neuron features to values between 0 and 1 . 

2. Find the similarity values between pairs of neurons using a 
predefined distance. 

3. Compute the preference values for each neuron. 

4. Cluster neurons using affinity propagation. 

5. Assign to all neurons the class determined by its exemplar. 

6. Compute the classification accuracy. 



To compute the similarity measure (step 2 in Algorithm 1), the 
Spearman distance, i.e., one minus the sample Spearman's rank 
correlation between observations (treated as sequences of values), 
was used. 

The similarity measure is computed as the opposite of the 
distance, s(i, k) = -d(i, k). Original settings of the affinity propa- 
gation algorithm are used (Frey and Dueck, 2007). The preference 
values (step 3 of Algorithm 1) for all points s(k, k) is computed as 
the median value of the similarity values s(i, j). 

To measure the agreement between a clustering produced by 
affinity propagation and a priori known set of labels for the 
points, a convention has to be set. We decide that the exem- 
plar of each cluster will determine the putative label of all points 
in the cluster. Consequently, the classification accuracy (step 6 
of Algorithm 1) is computed as the proportion of points whose 
putative label agrees with its true class. 
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We include a number of remarks concerning the application of 
Algorithm 1: 

- The output of affinity propagation (step 4 of Algorithm 1) 
depends on its input parameters. In particular it is sensitive 
to the similarity values between the points and the preference 
of each point to become a cluster. This means that changing 
these input parameters may determine changes in the number 
of clusters and their composition. 

- To evaluate the quality of a clustering produced by affinity 
propagation, two aspects should be simultaneously considered: 
(1) the number of points that are correctly classified. Since 
we assume that the labels of the exemplars are known, the 
correctly classified points may artificially increase. Therefore, 
we compute the classification accuracy as the ratio between 
correctly classified points (excluding the exemplars) and the 
total number of points (excluding the exemplars); and (2) the 
number of clusters, which should be preferably few. 

RESULTS 

DATABASE OF FOUR KNOWN INTERNEURON SUBTYPES 

We explored the use of affinity propagation to classify neo- 
cortical GABAergic neurons based on their morphological and 



Table 5 | Affinity propagation vs. Ward's method performance. 


Morphology Physiology 


Morphology and 
Physiology 


N clusters Accuracy N clusters Accuracy 


N clusters Accuracy 




AFFINITY PROPAGATION 


10 0.7374 36 0.8505 


8 0.7857 




4 0.5714 4 0.7575 
10 0.5859 36 0.8510 


4 0.6304 
8 0.6667 



W clusters is the number of clusters. Accuracy is calculated with respect to 4 
classes of neurons (BC, ChC, MC, non-MC). 



physiological properties. In order to test the affinity propaga- 
tion algorithm we used a well-characterized database where the 
group identity of the neurons was known from previous studies 
(McGarry et al., 2010; Packer and Yuste, 2011; Woodruff et al, 
2011). Specifically, we used a physiology database which con- 
tained 337 interneurons neurons distributed as: 63 BC, 218 ChC, 
40 MC, and 16 somatostatin positive interneurons non-MC. The 
morphology database contained 109 neurons distributed as: 31 
BC, 23 ChC, 33 MC, and 22 non-MC. Finally, there were 50 
neurons in a joint Morphology+Physiology database, formed by 
intersecting both databases. Its distribution was: 1 1 BC, 23 ChC, 
9 MC, and 7 non-MC. 

AFFINITY PROPAGATION CLASSIFICATION OF INTERNEURON 
MORPHOLOGIES 

Using these databases as ground truth, we performed affinity 
propagation, blind to the cell types. We analyzed the clusters 
structure, starting with the morphological database (Figure 1). In 
these figures, each neuron is shown as a colored glyph, i.e., a star 
plot that represents each neuron as a "star" whose ith spoke is pro- 
portional in length to its ith feature value. Similar shaped glyphs 
correspond to neurons with similar values of their features. The 



Table 6 | Affinity propagation vs. Ward's method performance. 

Morphology Physiology Morphology and 

Physiology 

N clusters Accuracy N clusters Accuracy N clusters Accuracy 



AFFINITY PROPAGATION 



10 0.8585 


36 


0.8471 


8 


0.9762 


WARD'S METHOD 


2 0.8037 


2 


0.9881 


2 


0.9792 


4 0.8000 


4 


0.9880 


4 


0.9783 


10 0.7879 


36 


0.9967 


8 


0.9762 



N clusters is the number of clusters. Accuracy is calculated with respect to 2 
classes of neurons (PV, SOM) 

















fi 
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FIGURE 4 | Hierarchical clustering found by Ward's method for the respectively represent neuron types BC, MC, non-MC and ChC. Neurons 

morphology and physiology database. Neurons are represented by are grouped into eight clusters and in each cluster the exemplar is 

their index in the database. Colors red, blue, green, and magenta, emphasized in bold. 
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colors provide additional information about the a priori known 
neuron class, i.e., red is BC, magenta is ChC, blue is MC and green 
is non-MC. Good clusters are those where neurons in the same 
cluster share the same color. BC and ChC (red and magenta) are 
closely related PV+ interneurons and MC and non-MC (blue and 
green) are both subtypes of SOM+ cells. 

The analysis of the morphology database resulted in 10 clus- 
ters (Figure 1). For example, the first cluster included 17 neurons, 
with 14 ChCs, had three misclassified neurons, those with num- 
bers 12 (of non-MC type, in green), 84 (MC, in blue), and 1 
(BC, in red). Its representative or exemplar was neuron 38 (a 
ChC). Clusters 3, 4, 5, 6, 7, 9, and 10 were also mixed, although 
in most of them, a single cell type also dominated. Clusters 2 
and 8 did not contain any error. To measure the accuracy of the 
classification, we counted the correctly classified neurons in each 
cluster. Overall, 73 neurons out of 99 are correctly classified (once 
the exemplars are not considered), yielding an accuracy of 0.73 
(Table 1). 

AFFINITY PROPAGATION CLASSIFICATION OF INTERNEURON 
PHYSIOLOGIES 

We performed a similar analysis of the physiological database, 
finding 36 clusters (Figure 2). Seven clusters (numbers 11, 17, 18, 
19, 20, 28, 33, and 35) are only composed by one cell. In con- 
trast, other clusters include many neurons. Moreover, some of 
these crowded clusters contained ChC neurons (magenta; clusters 
2, 9, 24, and 26), suggesting a potential method of identifying new 
subgroups of ChC cells, by analyzing some representative electro- 
physiological features (see below). Overall, the accuracy in this 
database was 0.85 (255 correctly classified out of 301; again with 
exemplars not considered in this calculation). 

AFFINITY PROPAGATION CLASSIFICATION OF INTERNEURON JOINT 
DATABASES 

We also performed an analysis of the combined anatomical and 
physiological databases, which had fewer neurons (50; Figure 3). 
We found 8 clusters, dominated by a single cell type. In fact, 
clusters 3, 4, 5, and 8 had no errors. From a total of 42 
non-exemplar neurons, 33 are correctly classified (accuracy of 
0.78). Remarkably, the simplest binary distinction parvalbu- 
min/somatostatin distinction was correctly identified in these 
clusters, with only 1 error out of 50 (neuron number 50, of BC 
type, in red). 

DEFINING CHARACTERISTICS OF THE 4 KNOWN INTERNEURON 
SUBTYPES 

As a final step, to extract further information from the clus- 
ters produced by the affinity propagation algorithm we identified 
those features that can serve to characterize each morphologi- 
cal and physiological cluster. The Wilcoxon rank sum test for 
equal medians was applied to compute features that have a sig- 
nificantly different distribution within each cluster with respect 
to the distribution in the whole database. 

Using this approach, a host of different variables were iden- 
tified, covering a wide spectrum of morphological and physio- 
logical features (Tables 2-4). These features could be useful in 
providing a compact characterization of those neurons included 
in the cluster. 



COMPARISON OF AFFINITY PROPAGATION AND WARD'S METHOD 

Ward's method of hierarchical cluster analysis is the current stan- 
dard used for classifying neuronal cell types. We tested whether 
affinity propagation improves accuracy over Ward's method. We 
performed the comparison at three different cluster numbers to 
illustrate the differences in performance. (1)4 clusters: To evalu- 
ate to what extent 4 clusters found by Ward's method correspond 
to the four known classes of neurons (PV-basket, PV-chandelier 
cell, SOM-Martinotti, SOM-non Martinotti). (2) The number of 
clusters found by affinity propagation (differs for each dataset): 
To evaluate the accuracy of Ward's method at the same number of 
clusters as affinity propagation. (3) 2 clusters: To evaluate to what 
extent Ward's method can separate neurons according to the 2 
largest neuron classes (PV vs. SOM). 

Accuracy was computed with the respect to the 4 class distinc- 
tions for scenarios 1 and 2 and with respect to the 2 class distinc- 
tion for all three scenarios. The exemplar for Ward's method is 
computed as the member of the cluster with the smallest mean 
distance to other members in its cluster. In every case the com- 
putation of accuracy was done without counting the exemplar as 
was done for affinity propagation. 

Table 5 shows the comparison of the two classification meth- 
ods, evaluated for 4 classes (scenarios 1 and 2). Figure 4 shows 
the hierarchical clustering found by Ward's method for the 
Morphology+Physiology database. Table 6 shows comparison of 
the two classification methods evaluated for 2 classes (scenarios 1, 
2, and 3). 

From these two tables we can see that affinity propagation per- 
forms better than Ward's method when accuracy is evaluated on 
4 classes, while the performance when accuracy is evaluated on 
2 classes is mostly similar between the two methods. Thus, affin- 
ity propagation can improve classification for finer distinctions 
in a dataset. This may be related to the capacity of affinity prop- 
agation to identify good descriptors of a subgroup of neurons 
represented by the exemplars and its ability to pass information 
between points in the dataset. 

DISCUSSION 

AFFINITY PROPAGATION: A NEW CLASSIFICATION METHOD FOR 
NEURONAL DATA 

In this methodological study we have explored the use of a new 
algorithm, affinity propagation, for the problem of quantita- 
tively classifying neuronal cell groups. We used a database of 
337 neocortical GABAergic interneurons, previously known to 
belong to four subtypes as "ground truth," and applied affinity 
propagation to a collection of morphological and physiological 
variables measured for each neuron. The classification accuracy 
we found is overall high: 0.73 for the Morphology database, 
0.85 for the Physiology database and 0.78 for the combined 
Morphology+Physiology database. The accuracy was higher 
when using physiological information from the neurons, but it 
contains a larger number of neurons than the other databases, 
so one would expect that the clustering algorithm will require 
a higher number of clusters to capture the diversity in the 
database. Moreover, grouping cells according to a binary parval- 
bumin/somatostatin class resulted in a much higher accuracy: 85 
out of 99 neurons were correctly classified in the morphological 
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database (accuracy of 0.86), and 255 out of 301 in the physio- 
logical database (accuracy of 0.85). In this respect, the perfor- 
mance of the affinity propagation algorithm in the combined 
Morphology+Physiology database was essentially perfect: 49 out 
of 50 neurons were correctly classified as being parvalbumin or 
SOM+. 

We conclude that affinity propagation can achieve a relatively 
high accuracy in the classification of interneurons, and is particu- 
larly accurate in distinguishing parvalbumin from somatostatin 
cells. It performs with higher accuracy than Ward's method in 
further distinguishing the four subtypes of parvalbumin and 
somatostatin interneurons. This unsupervised method, which 
already results in high classification accuracies, could improve 
if combined with dimensionality reduction, as we have found 
with other unsupervised or supervised neuronal classifications 
(Guerraetal, 2011). 

The extent to which it is useful to subdivide neuronal classes 
remains a debated subject. Some researchers view neuronal diver- 
sity as a continuum of features and refute the existence of discrete 
subtypes. However, affinity propagation can still be a useful 
method whether one is a "lumper" or "splitter" of neurons. 
Affinity propagation does not require an assumption that there 
are discrete subtypes for clustering. The clusters are produced 
using exclusively information about the features, not the num- 
ber of classes. We use the assumption of discrete numbers of 
neuron types to evaluate the output of the algorithms, but the 
results — relations between neurons — can be interpreted without 
this assumption. Additionally even if one takes the view of a 
continuum of neurons, the exemplars produced by affinity propa- 
gation are still valuable since they identify particular neurons that 
condense information about variations along the continuum. 

The identification of exemplars as a part of the clustering 
algorithm represents a major advantage of affinity propagation. 
Affinity propagation chooses an exemplar neuron for each cluster 
that can be studied as a representative of other neurons. The con- 
ventional construction of an exemplar for other clustering meth- 
ods is to compute the centroid of each cluster, computed as having 
the cluster average for each feature. Clearly these constructed 
neurons might not be realistic, while the affinity propagation 
exemplar is a real neuron from the dataset. 

While it is possible to pick a real neuron in each cluster as the 
exemplar as we did for our Ward's method analysis by choosing 
the neuron with the lowest mean distance to all other neurons 
in its cluster, affinity propagation has a more robust and infor- 
mative computation of the exemplar. By including the search for 
the exemplar as part of the clustering process affinity propagation 
does not restrict the set of potential candidate exemplars. When 
calculating exemplars after the clustering is complete, for Ward's 
method or any other method, the search for candidate exemplars 
and thus the candidate exemplars of members in that cluster, is 
restricted to the cluster. Affinity propagation takes into account 
the whole data set. 
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